#%%
import math
import os
import datetime
import pydicom
import matplotlib.pyplot as plt
import numpy as np

def gettimeformat(timevalue):
    return '%H%M%S.%f' if timevalue.find('.') > 0 else '%H%M%S'

def getdatetime(datestr, timestr):
    date = datetime.datetime.strptime(datestr, '%Y%m%d')
    time = datetime.datetime.strptime(timestr, gettimeformat(timestr))
    return datetime.datetime.combine(date, time.time())
    
def calculatedecay(actualframeduration, framereferencetime, halflife, acquisitiondatetime):
    '''
    定义计算衰变开始时间

    Parameters
    ----------
    actualframeduration : 
        (0018,1242) IS Actual Frame Duration
    framereferencetime : 
        (0054,1300) DS  Frame Reference Time
    halflife : 
        (0018,1075) DS  Radionuclide Half Life
    acquisitiondatetime :
        (0008,0022) DA Acquisition Date
        (0008,0032) TM Acquisition Time

    Returns
    -------
    scandatetime
        计算后的衰变开始时间，约为扫描开始时间
    '''
    frameduration =	actualframeduration / 1000
    decayconstant = math.log(2) / halflife
    decayduringframe = decayconstant * frameduration
    averagecountratetimewithinframe = 1 / decayconstant * math.log(decayduringframe / (1-math.exp(-decayduringframe)))
    scandatetime = acquisitiondatetime - datetime.timedelta(seconds = framereferencetime / 1000) + datetime.timedelta(seconds = averagecountratetimewithinframe)
    return scandatetime

#%% [markdown]
# ** 取第一张图进行测试数据，进行计算验证：**
# --- 
#| Tag        | VR   |  Value | Description |
#| --------   | :-----:  | :----  |  :----  |
#| (0018,1242) | IS | 150000 | # 1, 6 Actual Frame Duration |
#| (0054,1300) | DS | 74901.335499637 | # 1, 16 Frame Reference Time |
#| (0018,1075) | DS | 6586.2 |  # 1, 6 Radionuclide Half Life |
#| (0008,0022) | DA | 20170114 | # 1, 8 Acquisition Date |
#| (0008,0032) | TM | 145019.000013 | # 1, 14 Acquisition Time |


# %%
actualframeduration = 150000
framereferencetime = 74901.335499637
halflife = 6586.2
acqdt = getdatetime("20170114", "145019.000013")
scandatetime = calculatedecay(actualframeduration, framereferencetime, halflife, acqdt)
print(scandatetime)
#%% [markdown]
# --- 
# ** 期望结果：**
# --- 
#| Tag        | VR   |  Value | Description |
#| --------   | :-----:  | :----  | :----  | 
#| (0008,0021) | DA | 20170114 | # 1, 8 Series Date |
#| (0008,0031) | TM | 145019.000000 | # 1, 14 Series Time |
#%% [markdown]
# ---
# ** 读取整个PET序列并生成对比：**
# %%
def getdecaydatetime(ds):
    actualframeduration = float(ds[0x0018,0x1242].value)
    framereferencetime = float(ds[0x0054, 0x1300].value)
    halflife = float(ds[0x0054, 0x0016][0][0x0018, 0x1075].value)
    acqdt = getdatetime(ds[0x0008, 0x0022].value, ds[0x0008, 0x0032].value)
    return calculatedecay(actualframeduration, framereferencetime, halflife, acqdt)

def getserialdatetime(ds):
    return getdatetime(ds[0x0008, 0x0021].value, ds[0x0008, 0x0031].value)

def getacqdatetime(ds):
    return getdatetime(ds[0x0008, 0x0022].value, ds[0x0008, 0x0032].value)

def getdatetimes(imagedir):
    datasets = tuple([pydicom.dcmread(os.path.join(imagedir, filename)) for filename in os.listdir(imagedir)])
    scandatetimes = tuple([getdecaydatetime(ds), getserialdatetime(ds), getacqdatetime(ds)] for ds in datasets)

    count = 0
    for s, o, a in scandatetimes:
        if count % 30 == 0:
            print("index:{:0>3d},decaydt:{}, seriesdt:{} ".format(count, s, o, a))
            #print("index:{:0>3d},decaydt:{}, seriesdt:{},acq:{}".format(count, s, o, a))
        count = count + 1
    return scandatetimes

imagedir = ".\\pt\\siemens"
scandatetimes = getdatetimes(imagedir)

#%% [markdown]
# ---
# 生成时间曲线：

#%%
def showchart(scandatetimes):
    decaytime = [dt[0] for dt in scandatetimes]
    seriestime = [dt[1] for dt in scandatetimes]
    acqtime = [dt[2] for dt in scandatetimes]
    im = range(len(decaytime))
    plt.xlabel('Date Time')
    plt.ylabel('Image Index')
    plt.title('Calculate Decay Time')
    plt.plot(seriestime, im, 'green', label='Series DT')
    plt.plot(decaytime, im, 'blue', label='Decay DT')
    plt.plot(acqtime, im, 'red', label='Acquisition DT')
    plt.legend()
    plt.show()

showchart(scandatetimes)
#%% [markdown]
# ---
# 从曲线可以看出， 根据公式计算的衰变起始时间(蓝色)与序列时间(绿色)，基本一致，重叠在一起。采集时间是按照床位时间递进的。

#%% [markdown]
# ---
# 下面是我们的数据在吉大一扫描模体的PET序列验证。

#%%
imagedir = ".\\pt\\IQModel"
scandatetimes = getdatetimes(imagedir)
showchart(scandatetimes)

#%% [markdown]
# ---
# 下面是NueWise桂林临床的PET序列验证。

#%%
imagedir = ".\\pt\\NeuWise_GuiLin"
scandatetimes = getdatetimes(imagedir)
showchart(scandatetimes)
